%Take draws from the parameter space. 
%These draws are held in BETAS which is size NPxNVxNDRAWS
%Then calculate the probability of each person's sequence of choices at
%each draw for that person. 
%These probabilities are held in PROBS which is size NPxNDRAWS

BETAS=randi(NGridPts,[NP,NV,NDRAWS]); %BETAS go from 1 to NGridPts
BETAS=(BETAS-1)./(NGridPts-1);   %Now BETAS go from zero to 1 inclusive
for r=1:NV
  BETAS(:,r,:)=COEF(r,1)+(COEF(r,2)-COEF(r,1)).*BETAS(:,r,:);  %Now BETAS go from lower limit to upper limit for each coefficient
end;

rrate = squeeze(BETAS(:,1,:));
Psi = zeros(NP,T+1,NDRAWS);

if discountfun == "Exponential"
    for ii=0:T
        Psi(:,(ii+1),:)=(1./(1+rrate)).^(ii);
    end
elseif discountfun == "HM"
    for ii=0:T
        Psi(:,(ii+1),:)= 1./(1+(rrate.*ii));
    end
elseif discountfun == "Harvey"
    for ii= 0:T
        Psi(:,(ii+1),:)=(1./(1+ii)).^(rrate);
    end
%%% Quasi Hyper would take more set up work...
% elseif discount.fun =="QuasiHyper"
%    Psi(:,1,:) = 1;
%    for ii=1:T
%       Psi(:,(ii+1),:)= QH .* ((1 ./(1+rrate)).^(ii));
%    end
else
   disp('Discount.fun not found') 
end
    
%expand to correct number of rows using the ID column
Psi = Psi(XMAT(:,1),:,:);

% load values of discounted variables
%FirstRun will change to zero in bootstrapping, so rows align with
%bootstrap sample
if (FirstRun==1 && subset==0)
    pricetime = csvread('price.time.csv', 1)./10;
    captime = csvread('CAP.time.csv', 1);
    swtime = csvread('SW.time.csv', 1);
    buffertime = csvread('buffer.time.csv', 1);
    jobstime = csvread('jobs.time.csv', 1);
    infratime = csvread('infra.time.csv', 1);
    wlifetime = csvread('wlife.time.csv', 1);
    qualitytime = csvread('quality.time.csv', 1);
elseif (FirstRun==1 && subset==1)
    pricetime = csvread('price.time1.csv', 1)./10;
    captime = csvread('CAP.time1.csv', 1);
    swtime = csvread('SW.time1.csv', 1);
    buffertime = csvread('buffer.time1.csv', 1);
    jobstime = csvread('jobs.time1.csv', 1);
    infratime = csvread('infra.time1.csv', 1);
    wlifetime = csvread('wlife.time1.csv', 1);
    qualitytime = csvread('quality.time1.csv', 1);
elseif (FirstRun==1 && subset==2)
    pricetime = csvread('price.time2.csv', 1)./10;
    captime = csvread('CAP.time2.csv', 1);
    swtime = csvread('SW.time2.csv', 1);
    buffertime = csvread('buffer.time2.csv', 1);
    jobstime = csvread('jobs.time2.csv', 1);
    infratime = csvread('infra.time2.csv', 1);
    wlifetime = csvread('wlife.time2.csv', 1);
    qualitytime = csvread('quality.time2.csv', 1);
elseif (FirstRun==1 && subset==3)
    pricetime = csvread('price.time3.csv', 1)./10;
    captime = csvread('CAP.time3.csv', 1);
    swtime = csvread('SW.time3.csv', 1);
    buffertime = csvread('buffer.time3.csv', 1);
    jobstime = csvread('jobs.time3.csv', 1);
    infratime = csvread('infra.time3.csv', 1);
    wlifetime = csvread('wlife.time3.csv', 1);
    qualitytime = csvread('quality.time3.csv', 1);    
elseif (FirstRun==1 && subset==4)
    pricetime = csvread('price.time4.csv', 1)./10;
    captime = csvread('CAP.time4.csv', 1);
    swtime = csvread('SW.time4.csv', 1);
    buffertime = csvread('buffer.time4.csv', 1);
    jobstime = csvread('jobs.time4.csv', 1);
    infratime = csvread('infra.time4.csv', 1);
    wlifetime = csvread('wlife.time4.csv', 1);
    qualitytime = csvread('quality.time4.csv', 1);  
end


dbenefits = zeros(NROWS, 7, NDRAWS);
dprice = zeros(NROWS, NDRAWS);
for ii=1:NDRAWS
        dprice(:,ii) = sum(squeeze(Psi(:,:,ii)) .* pricetime, 2);
        dbenefits(:,1,ii) = sum(squeeze(Psi(:,:,ii)) .* captime, 2);
        dbenefits(:,2,ii) = sum(squeeze(Psi(:,:,ii)) .* swtime, 2);
        dbenefits(:,3,ii) = sum(squeeze(Psi(:,:,ii)) .* buffertime, 2);
        dbenefits(:,5,ii) = sum(squeeze(Psi(:,:,ii)) .* jobstime, 2);
        dbenefits(:,4,ii) = sum(squeeze(Psi(:,:,ii)) .* infratime, 2);
        dbenefits(:,7,ii) = sum(squeeze(Psi(:,:,ii)) .* wlifetime, 2);
        dbenefits(:,6,ii) = sum(squeeze(Psi(:,:,ii)) .* qualitytime, 2);
end

%Add up the non-price variables times their WTP
v=sum(bsxfun(@times,dbenefits,BETAS(XMAT(:,1),3:end,:)),2);      %sum(NROWSx(NV-1)xNDRAWS,2) gives NROWSx1xNDRAWS
v=squeeze(v);                              %NROWSxNDRAWS
v=bsxfun(@minus,v,dprice); %Subtract price
v=squeeze(BETAS(XMAT(:,1),2,:)).*v; %Multiply everything by price/scale coef
v=exp(v);
v(isinf(v))= 1.8e307;
sparsematrix=bsxfun(@eq,sparse(1:NCS)',XMAT(:,2)'); %NCSxNROWS
denom=double(sparsematrix)*v;                           %NCSxNDRAWS
denom(denom==0)=0.0000000000000001;
p=v(XMAT(:,3)==1,:)./denom;                     %NCSxNDRAWS
p(p==0)=0.0000000000000001;
personid=XMAT(XMAT(:,3)==1,1);                  %NCSx1
sparsematrix=bsxfun(@eq,sparse(1:NP)',personid');    %NPxNCS
PROBS=double(sparsematrix)*log(p);                      %NPxNDRAWS
PROBS=exp(PROBS);                               %NPxNDRAWS
clear sparsematrix v denom p personid dprice dbenefits Psi rrate
if WantBoot==0
    clear qualitytime wlifetime infratime jobstime buffertime swtime captime
end



